Gravitational radiation from radial infall of a particle into a Schwarzschild black hole. 
A numerical study of the spectra, quasi-normal modes and power-law tails. 
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The computation of the gravitational radiation emitted by a particle falling into a Schwarzschild 
black hole is a classic problem studied already in the 1970s. Here we present a detailed numerical 
analysis of the case of radial infall starting at infinity with no initial velocity. We compute the 
radiated waveforms, spectra and energies for multipoles up to Z = 6, improving significantly on 
the numerical accuracy of existing results. This is done by integrating the Zerilli equation in the 
frequency domain using the Green's function method. The resulting wave exhibits a "ring-down" 
phase whose dominant contribution is a superposition of the quasi-normal modes of the black hole. 
The numerical accuracy allows us to recover the frequencies of these modes through a fit of that part 
of the wave. Comparing with direct computations of the quasi-normal modes we reach a ~ 10~^ to 
~ 10~^ accuracy for the first two overtones of each multipole. Our numerical accuracy also allows 
us to display the power-law tail that the wave develops after the ring-down has been exponentially 
cut-off. The amplitude of this contribution is ~ 10^ to ~ 10^ times smaller than the typical scale of 
the wave. 



I. INTRODUCTION 

The gravitational radiation due to a point-like parti- 
cle falling into a black hole (BH) is a classic problem 
in General Relativity whose first computation goes back 
to the 1970s [U [5] . Being an elementary BH perturba- 
tion, it served as one of the first numerical applications 
of the BH perturbation theory initially developed in the 
classic papers [31 S] . The present literature in this field 
is abundanlj^ and the infalling particle model sometimes 
appears as a limit case of more general scenarios such as 
the infall of an arbitrary number of particles or, most 
importantly, the coalescence of BH binaries [5] , where it 
corresponds to the extreme mass-ratio limit. 

The general features of the radiated waveform are by 
now well understood [5]. The first part of the signal the 
observer receives is called the precursor and corresponds 
to the radiation emitted directly from the infalling source 
to the observer. It is therefore insensitive to what hap- 
pens near the BH and it was shown that it can be well de- 
scribed by a resummation of the Post-Newtonian expan- 
sion jlOHlSj ^ Then comes the ring- down phase which is 
dominated by a superposition of the quasi-normal modes 
(QNMs) of the BH. These are characteristic information 
of the Schwarzschild metric and are therefore brought by 
waves that were reflected in the neighborhood of the max- 
imum of the effective potential of the problem (the "bar- 
rier"), situated at ~ 1.5 Schwarzschild radii. The QNMs 
are a very interesting feature also because they provide 
a possible bridge between classical and quantum gravity 
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^ See [Sj for a summary of the theory and related references in 
theoretical and numerical studies. The textbook [6] also contains 
discussions and references on the subject. 

^ This was actually shown in the more general case of the coales- 
cence of a BH binary where the corresponding part of the signal 
is the oscillation due to the inspiral phase of the BHs. 



[Tll417j . Finally, the wave exhibits a power-law tail at 
large values of the retarded time, after the ring-down has 
been exponentially cut-off. This residual radiation corre- 
sponds to waves that were not initially heading towards 
the observer but ended up reaching him by scattering off 
the background metric at large distance of the horizon, 
hence the delay and the decrease in amplitude. 

In this paper, we focus on the simplest case, the radial 
infall starting at infinity with no initial velocity, as in the 
original works [H [2] (computations for more general ini- 
tial data include for example [3 [HI [19] ) . The spectrum 
of the radiated signal can be computed from a numerical 
integration of a single wave-equation, the Zerilli equation 
[3] . The first part of our work consists in performing this 
computation for multipoles up to Z = 6 with high ac- 
curacy, that is ~ 10~^. This allows us to analyze the 
spectra at quantitatively higher orders than present es- 
timates. The waveforms, which are obtained by Fourier 
transforming the spectra, reach an accuracy of ~ 10"'^ 
and ^ 10^^ at worst. For Z < 4, we are able to keep that 
precision on a relatively large interval of retarded time 
which includes the beginning of the power-law phase. 
Finally, we extract the QNM frequencies by fitting the 
ring-down phase with damped sinuses and then compare 
the results with the values obtained through direct com- 
putations [20l - f22] . This way we determine how "visible" 
the QNMs are. 

The organization of the paper is as follows. In section 
[njwe recall the main formulas describing the production 
of gravitational radiation from a radially in-falling test 
mass. In section HI we present and discuss our results. 



Finally, an account of computational details and error 
estimation is presented in the appendix. 



II. THEORETICAL BACKGROUND 

In what follows, m and M are the masses of the particle 
and BH, respectively. We use units G = c = 1 and the 
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following definition of the Fourier transform 



*,(r,a;) 
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'It: 



(1) 



The clear distinction in mass scales between the BH 
and the particle allows the problem to be treated within 
BH perturbation theory to first order (linearized Ein- 
stein equations). Therefore, the particle is a test mass, 
i.e. moving along the geodesies of the unperturbed 
Schwarzschild space-time, and produces a perturbation 
/ip^ (of the Schwarzschild metric) which is a test tensor 
field (on the Schwarzschild geometry) since we only keep 
first order terms of it in the Einstein equations. These 
can be reduced to two independent one-dimensional 
scalar wave-equations for each I multipole (in the ten- 
sor spherical harmonics decomposition), where I repre- 
sents the total angular momentum number. The wave- 
equation's effective potential term has no m dependence 
since the background metric is spherically symmetric. In 
the case of the radial in-fall of a particle, nor does the 
source term because of cylindrical symmetry. Thus the 
m ^ modes are not excited. At each I, these equations 
describe the dynamics of the two eigenstates of the parity 
operator, which don't mix because of the symmetry of the 
background metric under parity, and together fully deter- 
mine the 2'-pole of /i^;^. In the case of the radial infall 
of a particle, only the so called polar modetj^are excited 
and we are left with one equation, the Zerilli equation. 
Denoting its solution by ipi , the radial dependence of the 
2'-pole of hfj,i, is ^ {l/r)ipi(t — r) in the radiation zone. 
The Zerilli equation on the frequency domain is ^ 

dl^^i + {u:^ -Vi{r))^i^Si{r,uj) (2) 

where as usual — r+2M ln(r/2M— 1), and the effective 
potential is 

2A^(A + 1) + 3g^A^ + I (gf )^ A + I (g^)^ 
(Ar + 3M)2 
with A = (/-!)(/ + 2)/2. The source term is 



(3) 



Ar + 3M 



2A 



2M *a;(Ar 



(4) 



where T{r) is determined by the geodesic of the particle 
in the Schwarzschild metric, 



2M 



■ loe 



2 / r \3/2 2 f ^ 
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^r/2M + l) (y^72M-l)"' 
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^ They are the ones picking (—1)' under parity, i.e. the "true" 
scalar, vector, etc, as opposed to the "pseudo" ones picking 



(—1)'+^, called axial modes. 



The differential equation ^ is solved with boundary con- 
ditions of purely ingoing waves at the Schwarzschild ra- 
dius and purely outgoing ones at infinity 



lim 4(;(r^,a;) y4;_i„(w)e 
lim *;(^*7'^) ^/,o«t(w)e*' 



(6) 
(7) 



These correspond to the fact that the source is always 
localized in space and therefore GWs can only be emit- 
ted towards the infinities. In order to compute the radi- 
ated amplitude of the w-mode j4/_o„t(a;), it is convenient 
to use the Green's function method. Let yj^{r^,io) and 
y'^ (r* , uj) be the solutions of the homogeneous equation 
of ([2]) 

dly^ + {cj'~mr))y^ 







with boundary conditions 

lim y'i^[r^,^Lj) 



r*— f — oo 



lim j/,+ (r*,a;) 



(8) 

(9) 
(10) 



so that they match ^ and ^ in the final solution. At 
— 00 the potential vanishes and y^ tends towards the 
analytical form. 



lim ?/; (r*,a;) = a;(a;)e*"'"* + /3/(w)e 



(11) 



Thus the reflection and transmission coefhcients of Vi for 
a monochromatic wave coming from plus infinity are 



Ri 



ai{uj) 



Ti 



1 



(12) 



which we also compute for reasons made clear in section 
|A 1[ The Wronskian may be written Wi{u}) = 2i(jj(3i{uj) 
and, in the radiation zone, one obtains 

'^i,outir,uj) = lim «'/(r,w) 



r— >oo 
oicjr* /"OO 



Si{f^,uj)yi {f^,uj)df^ (13) 



2iu}Pi{uj) 

which by definition of the outgoing w-mode ([t]) gives 
1 



A. 



Lout 



(^) 



2iujj3i{(jj) J_ 



Si{f^,ui)yi {f^,uj)df^ . (14) 



Given the choice of normalization in ([I]), the radiated 
energy spectrum of the Z-mode is [1] [4] 



dEi 1 {1 + 2)1 , ,12 

da7=32^(^3^" "^'-*^")" ■ 



(15) 



The waveform is found using the inverse Fourier trans- 
form. Introducing the retarded time u = t — r* , we have 



1 



I, out 



{r,uj)e- 



2tt J-c 

OO 





^dw 



(16) 
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The last equality comes from ([7]) and the fact that tpi is 
real, i.e. A/^o„t(— w) = Ai^out{i^- So the procedure con- 
sists in computing y£ through eq. fe| with initial con- 



dition given by eq. ([9]) , and then evaluating ( 14 ) by ex- 
tracting /?/ {(jj) out of eq. and performing the integral. 
Once we have Ai^outi^), the energy spectrum is given by 



( 15 1 and the wave- function is computed through eq. ( 16 1 



III. RESULTS 

In this section we present the results of our numerical 
integration. A detailed account of the numerical proce- 
dure and error estimation is given in appendices |^ and 
[B} respectively. From now on we simplify the notation by 
dropping the "out" subscript in o«t(w) and ipi,outi^) 
and writing = dEi/duj. The figures are collected at 
the end of the paper. 



A. Analysis of the frequency spectrum 

The top and middle panels of figure [T] give the energy 
spectrum in both linear and logarithmic scales, the 
modulus and phase (of The shape 

of the energy spectra is in agreement with the result of 
ref. [I]. In table |!] we list the radiated energies for every 
2'-pole, that is 



E, 



!i{uj)duj . 



(17) 



To better characterize the energy spectrum, we have also 
computed the values (w*,/;*) at which /;(w) is maximal 
and the following quantities: 



El Jo 



u>fi{uj)duj , 



(18) 
(19) 



We also use ft{oj) = J2i<6 fii^) estimation of the 

total spectrum. 



~ 0.0104to^/M given in [T]. To a first approximation. 
El seems to follow the exponential trend proposed in ^] , 
which is ^ e~^'. We find that the form 



TABLE I: From left to right: Muhipole total radiated 
energy, energy spectrum's maximal value, corresponding 
frequency, mean value and standard deviation (for /;(w) 
seen as a distribution). 



El ~ 0.56(3);°-9(2)g-i.75(5)/^2/^^ 



(20) 



actually provides a better fit. Assuming that energies 
for ^ > 6 follow this empirical law, we get the order of 
the contribution of the neglected multipoles in the total 
energy, that is i?;>6 ^ 10~^m^/M, which is less than 
our numerical precision for Ef. As for the maximum of 
ft{oo), we find it is situated at Mu!^ = 0.3118(1) whereas 
[I] finds ~ 0.32. 
We now study the asymptotic behavior of the spectra. 



1. Low frequency limit 

In this limit we find that our numerical results are very 
well fitted by 

log\U^^ 1) - 0i(c^ = 0)1 ~af + 6flog(a;) (21) 



and 



log(/i(^«l))~af + bl\og{oj) 



(22) 



that is, 0/ and // have a power-law behavior. Within 
our numerical precision, we find that, for the multipoles 
Z = 2, . . . , 6 that we have studied, 4'i{u! — 0) is very well 
reproduced b}{^ 



= 0) 



7r(? -3) 
6 



(23) 



In order to obtain the coefiicients af ,bf ^a^ , bf with great 
accuracy it is necessary to perform the fit in the very low 
frequency region. In Table [Til we show these coefficients, 
obtained by performing a fit at 2Mw ^ 10~^, and we 
find that 



fiioj < 1) - w 
reproduces our results very well. 



21/3 



(24) 



I 






2Muj*i 


2M{uj)i 


2Mai (cj) 


2 
3 
4 
5 
6 


9.1368(9)10"^ 
1.1004(1)10"^ 
1.4947(1)10""* 
2.1380(2)10"'^ 
3.1423(3)10"'^ 


3.5943(4)10"^ 
3.3977(3)10"^ 
4.0757(4)10"'' 
5.3971(5)10""' 
7.5273(8)10"'' 


0.61992(9) 

1.0534(1) 

1.4685(2) 

1.8688(2) 

2.2726(2) 


0.5224(1) 

0.8747(2) 

1.226(2) 

1.582(3) 

1.941(4) 


0.1961(6) 

0.271(1) 

0.329(2) 

0.375(3) 

0.413(4) 


t 


1.0411(1)10"^ 


3.7744(4)10"^ 


0.6236(2) 


0.5723(1) 


0.2529(1) 



I 


4> 

«r 


bt 




bl 


2 


1.38(1) 


0.706(3) 


-1.84(6) 


1.321(8) 


3 


1.54(2) 


0.707(2) 


-4.40(5) 


1.989(7) 


4 


1.64(1) 


0.707(1) 


-6.92(4) 


2.653(6) 



TABLE II: The first few fitting coefficients in eqs. (|21) 
and (1221). 



For I = 2 we are even able to go down to 2Maj ~ 10 ^ 
and we get 



fi^2{i^ < 1) ~ 0.176(3) (2Mw)i =^32(2)^2 



(25) 



The estimation Et of the total radiated energy (bot- * Not directly computable because of the ui ^ term in the source 

tom left of table |l| is to be compared with the value t*^™' ^ut through extrapolation. 
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For each given I, the low frequency asymptotic behavior 
for fi(u!) and = 0) can be computed analyticahy be- 
cause it is due to the motion of the particle far from the 
horizon, where the trajectory can be well approximated 
by the Newtonian one and the GW emission can be com- 
puted using the multipole expansion in flat space. For 
I = 2 the contribution comes from the mass quadrupole 
and the computation is performed in ref. |23] (see also 
section 4.3.1 of ref. 61). The result is 



/;=2(a;< 1) 



2y/^ r^(l/3) 
3/ Svr 

0.1774 {2Mu} f^^ 



(26) 



so our numerical result ( 25 ) reproduces very well the ex- 
act analytic behavior. This is a significant check of our 
numerical procedure. We have also performed this an- 
alytical computation for I = 3 and again got agreement 
with eqs. (23) and (24). Both equations can be combined 
into 



Ai(c^«l)~a,(i^)('-3)/3^ 



(27) 



where a; is a positive real number. Finally, the 4>i{uj) 
variation, which is not computable analytically, is found 
to be constant with respect to / within our error margins 



0,(w<l)-(/.,(tj = O)-cj"-^°^(5) 



2. High frequency limit 



(28) 



For this limit the full relativistic treatment is necessary 
and there is no simple analytical expression to compare 
with. Here we only focus on /;. The top right panel of 
figure [l] clearly suggests an exponential cutoff. We find 
that the fitting form 



log /;(w > 1) ~ a; + biu + ci log(w) , 



(29) 



gives the best fit. Table [III] lists the resulting values for 
the parameters. 





ai 


hi 


Ci 


2 


5.34(2) 


-12.34(4) 


-3.03(8) 


3 


10.06(1) 


-12.26(3) 


-4.20(9) 


4 


14.81(1) 


-12.18(3) 


-5.42(3) 


5 


19.73(1) 


-12.17(4) 


-6.60(8) 


6 


24.81(2) 


-12.15(5) 


-7.7(2) 



TABLE III: Fitting coefficients in eq. (29) 



B. Waveforms and quasi- normal modes 

In figures [2] to [6] we show, for I = 2, . . . , 6, the am- 
plitude spectrum Ai{u}), the waveform ipi{u), the fit of 
its ringdown phase using QNMs and its power-law tail 
(compared to the obtained QNM fit). 



In the waveform plots, we observe that the number of 
significant oscillations in the ring-down increases with I 
while its typical length appears to be the same for all 
I. The precursor's typical length decreases with I. The 
typical amplitude decreases quite fast, as suggested by 
the empirical law in the corresponding energy ( 20 ) . 



As for the recovery of the QNM frequencies by fitting 
the ring-down phase, the procedure has already worked 
very well in numerical studies on the coalescence of BH 
binaries [131 [H] . It is therefore expected to work well in 
the case of the infalling particle model since it is a spe- 
cial case of the latter and numerically much simpler to 
treat. Table HV] lists the results of the fit for the first three 
overtones while table |V] lists the values obtained through 
direct numerical computations [20^22] (see appendix \A\ 
eq. ( |A3[ ) for the fitting form) . On the top panel of figure 
[7]we see that our error margins cover the expected values 
only for the first overtone (n = 1). Looking at the plots 
with the fitting curves, we see that the first mode's dom- 
ination interval is pretty much the same for all L On the 
other hand, the power-law contribution sets in later or is 
relatively weaker with increasing I. The high-Z graphics 
are indeed the ones were the QNMs fitted the wave best, 
so the QNM contribution becomes more "visible" with 
increasing I. Inversely, in the case of Z = 2, the power- 
law tail contribution sets in so early that it prevents us 
from even having a visible superposition of fit and data. 

Finally, concerning the tail of the wave, analytical 
studies show that it is actually a superposition of 
power-laws and that one has to go quite far in retarded 
time in order to see the leading term dominate. How- 
ever, at such high values our precision breaks down so 
we cannot perform any useful fit but we still display the 
graphic result for Z = 2, 3, 4 at lower u. The error on the 
1 — 5,6 cases is already notable at relatively low u. 



I 


n = 1 


n = 2 


n = 3 


2 


0.747(9) 


-i0.17(2) 


0.45(3) 


- j0.51(9) 






3 


1.198(2) 


-i0.18(l) 


1.104(8) 


- j0.54(l) 


0.93(4) 


-il.5(3) 


4 


1.616(1) - 


i0.186(6) 


1.580(8) - 


i0.552(6) 


1.34(3) - 


il. 140(4) 


5 


2.025(1) - 


i0.190(4) 


1.97(1) - 


i0.556(9) 


1.82(4) - 


i0.895(8) 


6 


2.424(1) - 


i0.190(4) 


2.378(8) - 


i0.572(9) 


2.22(4) - 


i0.805(4) 



TABLE IV: Computed QNM frequencies through 
waveform's fit. 



I 


n = 1 


n = 2 


n = 3 


2 


0.7473 


-i0.1779 


0.6934 - 


iO.5478 






3 


1.1989 


- i0.1854 


1.1653 - 


iO.5626 


1.1034 


- iO.9582 


4 


1.6184 


-i0.1883 


1.5933 - 


jO.5687 


1.5454 


- iO.9598 


5 


2.0246 


-i0.1897 


2.0044 - 


j0.5716 


1.9654 


- i0.9607 


6 


2.4240 


-i0.1905 


2.4071 - 


iO.5733 


2.3741 


- i0.9611 



TABLE V: Direct computation QNM frequencies. 
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1. The I3i{ij) parameter 



Consider equation ( 11 ). The convergence towards that 



asymptotic form being too slow, we use the next order 
terms (as in [l8]): 

2/j-(r, » 1,0;) ~ ai{io)e'^-' + pi{oj)e-'^^' + 

— {ji{u:)e'^''' + Si{oj)e-'^'-') (Al) 

We take 47 equally spaced sample points over 5 typical 
period^ 27r/a; on and plug them in (Al). Having 



four complex unknowns, this gives us an overdetermined 
linear system Ax = b, where A is a 47 x 4 complex matrix 
and X = {ai, Pij-fi^Si). Then multiplying by on the 
left we get a 4 x 4 matrix equation 

(A2) 



A^Ax 



A^b, 



the solution of which minimizes the residue \\Ax — 
This equation is then solved using the Gaussian elimina- 
tion method. 

This routine is used repeatedly at increasing , along 
with the computation of the integral in ( |14| , thus making 
the result more and more precisejj Note that there is no 
direct check on the error of /3i{uj). However, the fluctu- 
ations of its value strongly affecting Ai{uj), we take the 
latter's good convergence as a guarantee for an accept- 
able error on f3i{uj). A bit of monitoring at various stages 
of the computation confirms that the convergence of the 
(3i{uj) value is a lot faster than that of the integral in ( [l4| . 
Another source of confidence is the value Ri{uj) + TJJj) 
(see eq. ( 12 )) whose fiuctuation around 1 also gives qual- 
itative information on that error. In practice, we find a 
standard deviation of ~ 10~^ on the w-grid. 



Appendix A: Computational details 

In this section we present all the algorithms and tech- 
niques involved in the calculations. From now on we 
use only dimensionless variables. Thus, in what follows, 
t^, Ai^out, fu El and ijji^out actually stand for 2Muj, 
n/2M, Ai^out/{Mm), fi/m'^, M/rn^Ei and ipi^out/m, re- 
spectively. Relative errors and margins are denoted using 
the symbol S whereas A is used to denote absolute errors 
or integration grid steps, depending on the context. 



2. The radiation amplitude Ai(Ltj) 

For simplicity we write A for Ai{uj) since we describe 
the computation at a given value of w and I. We use Nu- 



These values are chosen so as to increase the information input 
in our fit and also avoid repeating values which lead to a badly 
conditioned system. 
® This is actually true only up to a certain limit because A 
contains three different orders of magnitude: (~ 1), (~ l/r*) and 
(~ so if r* is too large the problem is badly conditioned. 

In our case, however, we didn't reach that point. 
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merov's method to integrate in ([8|. Starting at finite 
r« <C 0, the initial condition corresponding to ^ is actu- 
ally e"*'^''* X]fc>o "fc(^^l)'° some constants = 0(1), 
at least for small k. We start at = —700 because that's 
where r — 1 starts being computable 10^"^°"), in dou- 
ble precision. This allows us to neglect the correction to 
the initial condition (|9|. Since the source term is also 
negligible at r* = —0(100), the integral in (14), which is 
computed in the same loop as i/j" , does not need any cor- 
rections for stating at finite either. We use the trape- 
zoidal rule for this integration because comparison with 
other Newton-Cotes formulae shows that it is the fastest 
method, in the sense that it starts approximating well at 
already big grid steps Ar* . This is due to the oscillatory 
nature of the integrand making positive and negative er- 
rors approximately compensate each other. The overall 
convergence is dramatically slow because the source term 
Q goes asymptotically as '-^ r~^/^. But the more we con- 
tinue the more the integration's error grows. Therefore, 
we have to use a special method in order to extract that 
limit faster. 



We define A{r^,) as being the value of (14) where the 
integral is truncated on the upper bound at r*, so that 
A = limr,_>.oo ^(t"*) (remember, the uj and I dependence 
is implicit here). On the actual r*fe grid (the discretized 
axis), the computed sequence Air^f^) approximating that 
value gives a damped oscillation. We divide the r*-axis 
into intervals of length L which is given by a few typi- 
cal periods 2tt/u! and index them by n = 0, 1, 2, . . .. We 
choose to consider only the average value of A(r*fc) out 
of every such interval, call it An (its phase </)„), a new 
sequence sharing the same limit. So L is a smoothing 
parameter. We also calculate /3; (and ai) at every n us- 
ing the procedure described in section |A 1[ The routine 
stops when the last ten values of |A„fand (/>„ are all 
within an e — 1Q~^ margin around their respective mean 
values and gives the latter as a result for A. The e mar- 
gin is actually a relative one for \A\'^ whereas it is an 
absolute one for cf). Thus, if we let the greatest distances 
from the mean values be denoted by A(|^p) and At/), the 
condition reads S{\A\'^) = A(|^p)/|Ap < e and A0 < e. 

Once that loop has finished, we start all over again but 
with half the previous Ar* step. This goes on until the 
difference between two consecutive such computations is 
again less than e (again, relative for \A\'^ and absolute for 
4>). When finished, we pass to the next w value. 

As for the w-grid parameters, let Aw/ be the step and 
i^i^max be the maximum value for which we perform the 
previous computation. It appears that uji^max ~ (^ + 1) 
is a good choice, as can be seen by looking at the Ai{Ld) 
plots (top left panels of figures [2] to |6]) where we have set 
{l + l)/2 for the maximum of the displayed uj axis. This is 
why we chose Aw/ — 6(^ + 1)10^^^, so that the precision is 
the same for all I. However, we did not choose the uJi^max 
value to follow the ~ (Z -I- 1) trend. We find instead that 
there is a natural limit on the w axis for the convergence 
of the Ai{uj) computation, given our precision criteria. 
After a given value for //, apparently common to all I 



the program keeps dividing the Ar* step without ever 
meeting the required precision. Since the values at a; 3> 
1 are important for the computation of the tail of the 
waveform, we set the LUi^max value the higher we can, 
that is to that natural limit. This also sets Z = 6 as our 
limit for /, because the computation for I = 7 would not 
give enough points for the spectrum at large frequencies. 



The radiated waveform tpi{u) and energy 
spectrum integrals Ei,{uj)i and (7;(cj) 



We consider equations ( 16 ) , ( 17 ) , ( 18 ) and ( 19 ) . We 



use a high order Newton-Cotes formula for their integra- 
tion, the one with 7 stages and of order 8, named Wed- 
dle's formula. We also use Richardson's extrapolation on 
the integrals computed with steps Aw and 2Aw in order 
to further increase our precision. 

For V'i(u), the / = 2 case must be treated carefully 
because Ai^2i^) diverges at the origin. It is true that, 
for any value of I, we can only compute Ai{uj) for lo > 
because of the w~ 



factor in the source term. However, 
we know from (27) that Ai=3{uj = 0) is finite and there- 



fore deducible by extrapolation and A/>3(a; = 0) = 0. 
For I = 2, the [0, Auj] contribution cannot be extrapo- 
lated. It can however be approximated through the inte- 
grand's analytic behavior at w — >■ (discussed in section 
III A neglecting the (j){u}) variation). Plugging the later 

we obtain a solution involv- 
ing generalized hypergeometric functions and we Taylor 
expand them until the desired precision is reached. 



4. Extraction of the QNMs out of the signal 

We know that the QNMs are getting more damped 
with increasing n. Thus, if we look far enough in the 
ringdown phase of the wave, the least damped mode (n = 
1) should dominate. However, the more we go at large u 
to look for the first mode the more the tail contribution 
becomes notable. So our fitting model is 



Si,«e'^''"-^" sin(i:j;^„^sRU -|- 6';,„) + /l;,„ 



(A3) 



where the last term is an offset which can help compen- 
sate the shift due to the slowly appearing tail. It is re- 
ally necessary for low I but becomes negligible for high I. 
Once we have found the first overtone (n — 1), we sub- 
tract it from the waveform and reapply the fit seeking 
the second one {n — 2) and so on. 

Since there is no prescription for finding the ideal in- 
terval on the u axis to perform the fit we run a small 
program which tries all possible intervals in [0,30], re- 
taining every time the value of the resulting fitting pa- 
rameters. This gives us a 3D plot for each one of them, 
where the floor axes are given by the values of the inter- 
val's boundaries. The plots for wsr and wq exhibit some 
flat areas at equal height corresponding to the ensemble 



7 



of intervals at which the fit was optimal. We then cut 
the plot in horizontal slices of a given thickness Az and 
create a histogram giving the number of points which lie 
inside each one of those slices. The maximum of that 
histogram gives us the researched value with absolute er- 
ror Az/2 but when the spike is not clear enough we take 
half its width for an estimation of the error instead. For 
n > 1 we sometimes obtain more than one maximum so 
we choose based on the coherence with the rest of the 
data and by looking at the 3D plot in order to identify 
"false" flat areas. The bottom of figure [7] gives an exam- 
ple of the 3D graphic and its histogram. As for the B 
and 6 parameters in ( A3 ) , no such flat area is obtained 



after the first fit series so we run the program one more 
time but with the wsr and ujcg values already inserted. 



Appendix B: Error estimation 
1. Radiated spectrum Ai{uj) and fi{u}) 



Almost all th e com puted points meet the precision cri- 
teria of section A 2 meaning an estimated ~ 10^^ pre- 
cision on \Ai\'^ and <j> and therefore on \Ai 
member that it is an absolute error for (h 



^ 10- 

and fi (re- 
The only 



and 

it is an absolute error 
exceptions are for / > 4 where points with w close to zero 
(a; ~ IQ-^) have an error of 10"^. However, \Ai\ is 
very small there and the absolute error it causes in the 
calculation of Ei and ipi is thus negligible. In order to 
decrease it as much as possible anyway, we have extrap- 
olated those regions using the analytical low frequency 
behavior discussed in section UlIAl 

As for the asymptotic behavior fits, in most cases such 
as in the high frequency region, there is no prescription 
for the ideal interval to fit on, so we estimate the error 
on the fitting parameters by their variation when fitting 
on different intervals. If, however, the number of points 
is small, leaving no choice about the fitting interval, the 
error corresponds to a 95% c.l. 



2. Radiated energy spectrum characteristics 
(table ^ 

There are three types of error in the computation of 
the integrals Ei , (w) ; and cri (uj) : the one due to the rela- 
tive error on the integrand, noted 6i, the one due to the 
discreteness of the integration domain, noted Sd, and the 
one due to its finiteness, noted Sf. Sd is estimated by 
the relative difference of two integrations with grid steps 
Aw and 2 Aw and Sf is estimated using (29). For Ei we 
have Si ~ 10"^ since ~ |Ai(a;)p, ~ 10"'^ and 

6f ^ 10"^ so S{Ei) ^ 10^^. For (a;); and o-;(w), we use 



Ap{xi, ...,Xn) 



k=l 



dp 



dxi. 



Axk 



for the Si 
2.10-4(1 - 



and find Si{{uj) 



2.10- 



The Sd and S 



'4 and St{ai{uj)) = 
are again smaller. 



Finally, for the maximum we record the smallest gap 
between its value and its direct neighbor's on the a; grid. 
However, this value is inside f*{l ± e), so the relative 
error for is also given by e. For the error on W;* we 
simply take half the uj grid step. 



3. Radiated waveform ipi{u) 

First of all, being interested in the domain u G 
[-50, 200] for ; < 4 and [-50, 50] for / = 5 and 6, the 
greatest gap between two consecutive lou values in the 
oscillatory term of ( 16 ) is {uAuj)max = 0.06 ^ tt. Thus 
the w-grid is dense enough to take into account even the 
sharper variations of the integrand. The sources of error 
are the same as in the previous section although here we 
are going to use the absolute analogues for Sf and Sd- The 
relative ones vary a lot near the zeros of the waveform 
and aren't therefore very meaningful. A/ is estimated 
using (29) and (15) to obtain the behavior of |Ai(ti;)| 



1)1 

Then 



'I287r 



(/ + 2)! 



-b,uj)/2^c,/2-l 



(B2) 



/27r 



\Ai{u! > 1)1 cos(uw)dw (B3) 



is a good estimation of the absolute error due to the 
neglected part of the frequency domain. Table |VI| shows 
the typical values for A^; and Af which fluctuate well 
inside the given order. It also gives the maximum values 
of the corresponding waveforms in order to compare the 
scales. 



/ 


Ad 


A/ 


|V'|i,max 


2 


10-** 


IQ-^ 


4.2 X 10-' 


3 


10-** 


10-9 


3.8 X 10-^ 


4 


io-« 


IQ-* 


6.3 X 10-^ 


5 


10-^ 


10-10 


1.2 X 10-^ 


6 


10-1° 


10-" 


2.6 X 10-* 



TABLE VI: Typical values of A^; and Af and reference 

scales 



The regions where the relative error is maximal are of 
course the ones near zeros and towards the end of the tail 
(~ 10-^ at worst for the latter). Otherwise, the majority 
of points has ^ IQ-^ or even 10^^. Finally, the erro r du e 
to the one of the integrand is Si ^ 10"^ (see section A 2 ). 

As for the QNM, we have already explained how the 



(Bl) errors are computed (see section A4) 
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FIG. 1: Top panels: the energy spectra l/m?dEi/duj on normal and logarithmic scale. Middle panels: The 
amplitude spectra, modulus (logarithmic scale) and phase. Bottom panels: low frequency asymptotic behavior of 

the energy spectra and phases. 



10 



CM 

3 




is- 






CM 


CM 


CM 


o 


OJ 


CM 


CM 


CM 


o 


O 


o 


O 


o 


o 


O 


O 


o 


1 


1 


1 


1 


+ 


1 




1 


1 


El] 






H 










W 




ro 


CM 




o 




Csl 


CO 















LO 


o 


IT) 




o 


o 


O 


o 


o 


o 


o 


o 


o 


1 


1 


1 


1 


1 


1 


+ 


1 


1 


El] 


ti] 


Ell 


PJ 


Ell 


EJ 


El] 


El] 


El] 


o 


LT) 


o 


LT) 


o 


o 


o 


o 


o 


ro 


CM 




t— 1 


iH 


LO 


o 


LO 


1 








Osl 


e 




tH 


M 


II 


II 









LH 


c 


G 








> 


2 


S 




Z 


z 




a 


a 



i^ 




CM 


CM 


CM 


CM 


o 


CM 


CM 


CM 


CM 


O 


O 


O 


O 


o 


O 


O 


O 


O 


1 

Ei] 


1 

ti] 


1 

Ei] 


1 

ti] 


+ 


1 

Ei] 


1 

ti] 


1 

Ei] 


1 

Ei] 




CO 


C\] 


t— 1 


o 




CM 


CO 





CD 

o .2 



as 

^ o 

CD U 
^ CD 

^ ^3 



(D 



a 

o 
o 



S o 

=2 CD 

S -a 
1 ^ 

!^ CD 

CD 

a 

,S O 

bb +^ 
'S o 

rO . 



«=; CD 
1 o 



CD 

-a 

pi 



1 

CD 

-13 



a 
o 



o 

I— I 

fa 



12 



a 

CD 



-5- 




ro 


00 




o 




00 


CO 


o 


o 


o 


o 


o 


o 


o 


1 


1 


1 


+ 


1 


1 


1 




Ci] 


t1 






Ci] 




LO 


o 


o 


o 


o 


o 


LO 


iH 


■—1 


LO 


o 


LO 


t— 1 







ro 


ro 




o 




CO 


00 


o 


o 


O 


o 


o 


o 


o 


Ci] 


1 

H 




+ 




1 




m 


o 


o 


o 


o 


o 


LO 


>— 1 




LO 


o 


LO 


iH 


1 



Pi 
o 



o 
o 

o 

CD 
,£3 




Pi 
ce 

a 



o 



o O 

-9- CD 

s § 

o 

O 
<D 
t/} 

(D 



0) 



H o 
. o 



CD 

pi 



CD 



its 

CD 
I — I 



O 
I— I 
fa 



13 




CZ) O CD O CZ) o o 5h 

ooooooo <D 

o o o o o o o ' ' 

ooooooo Q 

III ^ 



d 

I— I 

fa 



14 



2MbJn.Q VS. 2Maj„,sR 



-0.1 



-0.3 - 



-0.5 - 



-0.7 



-0.9- 



-1 . 1 



-1.3 



-1.5 



— --1 


= 2 


1 


= 3 


1 


= z 




1 


= 5 


---1 


= 6 



0.0 0.5 1.0 

Bi^n/Mm vs. overtone number n 



1.5 2.0 2.5 3.0 

Oi^n VS. overtone number n 



1 = 2 

1 = 3 

1 = 4 

1 = 5 

1 = 6 



-1 = 2 

1 = 3 

1 = 4 

■1 = 5 

■1 = 6 



tiJ3,i,K VS. [a, h] 



|{ points in the fc-th slice }| vs. k 




FIG. 7: Top panel: the computed QNMs through the waveforms fit. The dotted hues' intersections are the expected 
values. Vertical lines link QNMs of the same / while horizontal lines link QNMs with the same n. n increases by 
going downwards. Middle panels: the amplitudes Bi^n/Mm and phases 9i,n of the first two overtones. Bottom left 
panel: example of the 3D plot obtained when fitting the waveform for QNMs on all intervals [a, b] in [0, 30] of u/2M. 
Bottom right panel: the corresponding histogram giving the number of points in each horizontal slice of the 3D plot. 



